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ABSTRACT 

We present a 190-307 GHz broadband spectrum obtained with Z-Spec of NGC 1068 with new 
measurements of molecular rotational transitions. After combining our measurements with those 
previously published and considering the specific geometry of this Seyfert 2 galaxy, we conduct a multi- 
species Bayesian likelihood analysis of the density, temperature, and relative molecular abundances 
of HCN, HNC, CS, and HCO+. We find that these molecules trace warm (T > 100 K) gas of H2 
number densities 10^'^ — 10^'^ cm~'^. Our models also place strong constraints on the column densities 
and relative abundances of these molecules, as well as on the total mass in the circumnuclear disk. 
Using the uniform calibration afforded by the broad Z-Spec bandpass, we compare our line ratios to 
X-ray dominated region (XDR) and photon-dominated region models. The majority of our line ratios 
are consistent with the XDR models at the densities indicated by the likelihood analysis, lending 
substantial support to the emerging interpretation that the energetics in the circumnuclear disk of 
NGC 1068 are dominated by accretion onto an active galactic nucleus. 

Subject headings: Galaxies: Individual: NGC 1068 - Galaxies: Seyfert - Galaxies: ISM - ISM: molecules 
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1. INTRODUCTION 

Molecular transitions in the millimeter and submil- 
limeter bands are well suited to characterize the phys- 
ical conditions of dense gas in active galactic nuclei 
(AGN) or starburst (SB) regions. Starburst regions are 
associated with far-UV emission from O and B stars 
and corresponding p hoton-dominated regions (PDRs, 
iTielens fc HoUenbachl (|l985. '). which characterize cloud 
surfaces) , while accretion onto an AGN cr e ates a n X-ray 
dominated region (XDR, IMalonev et all (|T996l ). which 
characterize cloud volumes). These two different types 
of regions can be distinguished by the signatures of 
the molecular gas chemistry; making such a distinction 
has been the focus of recent research in many (ultra- 
)luminous infrared galaxies or (U)LIRGs. However, we 
often must understand the physical environment of the 
dense gas before we can properly diagnose the XDR or 
PDR nature of it. For example, the predicted line ratios 
for XDR or PDR scenarios vary with density. Therefore 
it is important to understand the physical conditions of 
the gas before using (sub) millimeter lines to provide key 
information about the dominant energy source heating 
the gas. 

NGC 1068 (M77) is often cited as the prototypical 
Seyfert 2 galaxy; it is a barred spiral hosting an AGN 
that has been extensively studi ed at many wavelen gths. 
At a redshift of 1137 kms'^ (|Huchra et al.llT999i) . 1" 
« 78 pc assuming Hq = 70 kms~^ Mpc~^. With an in- 
frared luminosity of 2 x 10^^ L0, it is classified as a LIRG 
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(jSanders et al.ll200l . The major features of NGC 1068 
are the large star-forming arms, a central circumnuclear 
disk (CND), and a jet emerging in the northeast direction 
from the CND. The star-forming arms, seen most promi- 
nently in CO line emission, essentially for m a large ring of 
30'^ di ameter (abo ut 2.3 kpc. Figure 2 of iHelfer fc Blitd 
(fT995l) . Figure 1 of ISchinnerer et al.l (12000|)). The CND 
is approximately 4" diameter (312 pc), from which arises 
emission of molecular high-density gas tracers (e.g., CS, 
H CO^, HCN, HN C ), as c an be seen in m any ma p s (suc h 
as IHelfer fc Blitz! (fl995l): iG arcia-BuriUo'^aLl (pOOl : 
IKrips et al.l (|2008[ ): lBavet et a.1., (.2009.) and the 1 mm 
and 3 mm continuum maps of IKrips' "etall (12001) 1. The 
CND itself can be resolved into two knots (East and 
West) with different velocity profiles (Uscro et al. 2004). 
The geometry of NGC 1068 is a significant considera- 
tion for single-dish observations, because the two very 
different regions (the CND and SF ring) are blended by 
most (sub)millimeter telescope beams. For example, for 
our ^ 30" beam, we are measuring some emission from 
the SF ring on the outer edges of our beam, and also 
emission from the central CND concentrated in the mid- 
dle. In the case of CO (and the continuum) we are mea- 
suring both components. For the high-density tracing 
molecules, however, the emission is concentrated only in 
the center of our beam (from the CND), as can be seen 
in the aforementioned interferometric maps. The physi- 
cal area that we are modeling for this paper is the CND 
only; we address how we correct for possible contamina- 
tion from the starburst ring in §3.41 

The black hole mass at the center of the CND, as 
inferred from water mase r emission, is 1.7 x 10'' M© 
(|Greenhill fc GwinnI I1997D . The bolometric luminosity 
of the entire galaxy (ring + CND, 2.5 - 3 x 10" Lq) 
is determined mainly by the mid-IR; about half of this 
emission comes from the starforming ring, so the bolo- 
metric luminosity of the cen tral 4" in diamete r has been 
estimated as 1.5 x 10" L^ (|Bock et al.|[2000l ). Though 



2 



KAMENETZKY ET AL. 



2004 : 



2008 



this is the luminosity visible to us from the CND, the 
AGN itself has a higher intrinsic luminosity that is well 
shielded by i ts obscuring torus. Using [OIV] 26 /xm as 
a calibrator, iRigbv et alJ (|2009f ) infer an intrinsic (but 
highly obscured) AGN luminosity of (3.1 ± 0.9) x 10" 
Lq, consistent with the estimation of iBock et al.l (|2OOO0 
(3 X 10^^ L0) based on Lboi ~ 30 lj2-wkeV- The high 
amount of the AGN's luminosity that is absorbed indi- 
cates that the AGN can be extremely important in pow- 
ering the CND emission. 

Several recent studies indicat e that the CND 
of NGC 1068 is an XD R (lUsero et all 
Perez-Beaupuits et al.l 120091 : iKrips et al.l 
Garcia-Burillo et al. r i2oion . There are a variety of 
diagnostics of XDR/PDR environments that can 
be utilized with millimeter-range molecular tran- 
sitions. AGN tend to have lower intensity ratios 
of HCO+/HCN (J = 1 0) and higher ratios of 
HCN/CO (J = 1 ^ 0, IKrips et al.l (gOOS)). HCO+ 
[(J = 3^2)/(J = 1^0)], HCN [(J = 2^1)/(J = 1^0)], 
and HCN [( J = 3 ^- 2)/(J = 1^-0)] intensity ratios also 
appear lower in AGN than in SB. In Krips' study of 12 
nearby galaxies with varying degrees of AGN/starburst 
activity, NGC 1068 was the most characteristic of AGN 
line ratios (and generally appearing completely opposite 
from M82, a prototypical starburst galaxy). Further- 
more, it is predicted that CN is more abundant than 
HCN in XDR scenarios due to dissociation (HCN 
H + CN, but then CN is parti c ularly r obust to further 
dissociation, iLepp fc Dalgarnol (|1996l ): iMeijerink et al.l 
(|2007)). 

By observing NGC 1068 with Z-Spec, a broadband 
millimeter- wave spectrometer, we are able to measure 
the continuum simultaneously with multiple molecular 
transitions. We present observations of HCN, HCO+, 
and HNC which complement the current literature as 
well as new transitions of CS. The paper is organized 
as follows. The Z-Spec instrument, our observing pro- 
cedure, and the spectrum are described in fj2] We use 
the measurements in our spectrum in conjunction with 
data from the literature to model first the physical con- 
ditions of the r uolecular gas in the CN D of NGC 1068 
using RADEX (|van der Tak et al.l 120071 ). The modehng 
procedure is described in ^21 with discussion of the re- 
sults beginning in 21 We next use the derived molecular 
abundances and the line ratios within Z-Spec's band to 
demonstrate the XDR nature of the CND in S|5l followed 
by conclusions in Sj6l 

2. OBSERVATIONS WITH Z-SPEC 

In this section, we first describe the Z-Spec instrument 
and observations of NGC 1068. We then present our 
spectrum along with a description of our spectral fitting 
procedure and a short discussion of the fit results. 

2.1. The Z-Spec Instrument 

Z-Spec is a broadband (190-307 GHz) millimeter- 
wave grating spectrometer which we have used at the 
Caltech Submillimeter Observatory (CSO). Its large 
bandwidth allows simultaneous observations of multiple 
molecular rotational transitions along with the under- 
lying continuum. The detector array is composed of 
160 silicon-nitride micromesh bolometers cooled to 60 



mK by an adiabatic demagnetization refrigerator (ADR) 
and a closed-cycle "^He refrigerator. Z-Spec's compact 
design is acheived via a WaFIRS (Waveguide Far IR 
Spectrometer) design utilizing a parallel-plate waveguide 
and curved Rowland diffraction grating (Bra dford et al.l 
120031) . Z-Spec's spectral resolution is approximately 900 
MHz at the band center, but varies from 500 MHz (700 
kms~^) at the low frequency end to 1300 MHz (1200 
kms~^) at the high frequency end of the band. To min- 
imize susceptibility to atmospheric fluctuations and to 
subtract sky background emission, we use a chop-and- 
nod mode. For NGC 1068 observations, our nod interval 
was 20 seconds, and the chop frequency was 0.95 Hz with 
a throw of 90". More de t ails a b out the instrum e nt ca n 
be found inlNavlor et all (l200l. 'Bradford et al.l (pOOl . 
lEarle et al.l (!2006( ). and llnami et a l. (200^7^ 

We observed NGC 1068 over three nights in late Jan- 
uary 2007 for a total of 4.21 hours of integration time 
(including chopping/nodding, meaning one half of this 
time was on source). On January 24 and 28, the optical 
depth at 225 GHz, T225GHZ was steadily ~ 0.05, while 
on January 27 it was slightly worse at T225GHZ ^ 0.070- 
0.086. The median sensitivity (channel uncertainty times 
square root of integration time) was 1.13 Jy s^/^, some 
two to three times greater than we typically achieve in 
observations of faint sources. We attribute the excess to 
systematic effects due to the bright source. 

Z-Spec's calibration relies on using a set of planetary 
observations with varying observing conditions (bolome- 
ter loading and bath temperature) and fitting the de- 
pendence of each individual bolo meter's respons e on it s 
operating voltage, as described in lBradford et al.l (|2009). 
Because our observations of NGC 1068 were conducted 
before the calibration curves described in the aforemen- 
tioned paper, it was more appropriate to calibrate these 
observations with local-time observations of Uranus. The 
total integrated continuum flux measured with the Brad- 
ford and Uranus calibrations agree to within 5%. 



2.2. Continuum Analysis and Line Fitting Procedure 

The Z-Spec spectrum of NGC 1068 is shown in Fig- 
ure [T] Our line-fitting procedure treats each emis- 
sion line as a single-component Gaussian. Each chan- 
nel's spectral response profile has been previously mea- 
sured and is incorporated into the line-fitting routine. 
Higher-resolution spectra reveal more detailed structure, 
such as two or more components for the high-density 
tracers and generally three c omponents f or CO (e.g. 
iPerez-Beaupuits et al.l (|2007f) : IKrips et al.l (|2008D ). but 
we cannot resolve these structures. Because Z-Spec can- 
not resolve line widths below ~ 1000 kms~^, the line 
width was set to 240 kms"""^ based on the approximate 
CO J = 2 ^ 1 linewidth in [Israel (l2009l) ( some mea- 
sured NGC 1068 linewidths have been higher or lower). 
One channel (at 213.3 GHz) was not operating during 
our observations and is excluded from our fit. In our 
line fitting procedure, we choose the lines to be fit based 
on a list of the expected transitions in our band, and 
those which are not reliably detected are excluded and 
then the spectrum refit. Rest line frequencies are taken 
from the online Leiden Atomic and Molecular Database 
^Schoier et al... 20051 
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Fi g. 1 . — Spectrum of NGC 1068 observed by Z-Spec on the CSO, shown in logarithmic (a) and linear (b) scale. The spectral fit (described 
in i|2.2l l is overplotted in green on the linear plot. Vertical fines indicate the frequencies of known molecular transitions of importance, 
including the prominent CO J = 2 — > 1 line (vertical lines alternate between red and blue for clarity only). Transitions of SiO are not 
included in the fit but are shown here for reference. The results for HOC+ J = 3 — > 2 and C'^^S J = 4 — > 3 are only upper limits. The 
unlabeled line at ~ 290 GHz is an unidentified feature that is modeled as a Gaussian in the fit. Features above CS J = 6 — ^ 5 are not 
included due to the high channel-to-channel variation at the end of t he b and due to atmospheric noise. The features in the spectral fit at 
200 and 260 GHz are instrumental artifacts described at the end of i|2.2l Panel (c) illustrates the difference (residual) between our model 
fit and the spectrum, with uncertainties. Panel (d) shows the contribution from each channel to the total x^; the ten channels around the 
CO J = 2— >1 line are responsible for 40% of the total of 5.1. 



iKrips et aLl (|2006[ ) measured the continuum emission 
from the core and jet of NGC 1068; at 230 GHz, the 
measured continuum flux from this region is only ~ 
10% of our measured continuum. Therefore, the con- 
tinuum measured by our beam is dominated by the 
starburst ring, but a small contribution is also present 
from the central disk and jet. We add their "core -I- 
jet" model (their Figure 3, top panel) to the (beam- 
sc a,led) 34 K greybody o f the starburst ring measured 
bv lSpinoglio et all (|2005[ ). so that our continuum fit be- 



comes: 



F,, =A 



240GHz 



B-2 



nB^{T)[ 1 - e 



~t~ -^O,core-t-jet 



V 230GHz 



(1) 



Jy- 



The first term is the beam-scaled greybody, and is the 
source size in steradians; we estimate fl to be approxi- 
mately 1.66 X 10~^ sr, based on the 30" diameter of the 
starburst ring. Since the emission from the ring actually 
only occupies a portion of this total area, this is likely an 
overestimation that can be accounted for with the free 
parameter A. B,j{T) is a blackbody of dust temperature 
T = 34 K, and we take vq = 3000 GHz (which corre- 
sponds to Ao — 100 /xm) and (3 ^ 2 fo r the dust emissivity 
spectral index (|Spinoglio et al.l[2"005[ ) . The second term is 
the power-law contribution from the disk and jet, where 
Fo,core+jet IS the mcasurcd flux at 230 GHz and a is the 
measured power-law scaling. We use J o.eore+ iet = 0.028 
Jy and a = 0.9 from iKrips et all (|2006D .' 

A and B are the free parameters in our fit; B is 
meant to represent the beam-scaling, because our beam- 
size changes with frequency and thus our coupling with 
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TABLE 1 
Z-Spec Detections in NGC1068 



Transition Urest f F„dv J TAdv S/N^ 

[GHz] [Jy kms-i] [K kms"!] 

C02^i 230.54 8366 183.7 436 

i3C02^i 220.4 712 15.8 54.4 

Ci^02^i 224.7 31 0.7 3.0 

Ci*02^i 219.6 236 5.2 16.6 

CCH3_>2 262.25 57 1.2 6.1 

CN2^i 226.87 496 10.9 54.2 

CS4^3 195.95 68 1.6 3.6 

CS5^4 244.94 58 1.2 5.4 

CSe^s 293.91 58 1.2 3.0 

C34S5^4 241.02 40 0.9 3.8 

HCN3^2 265.89 406 8.5 33.1 

HCO+3^2 267.56 267 5.6 15.9 

HNC3^2 271.98 104 2.2 10.3 

Unidentified'' 290.02 93 1.9 6.0 

C34S4^3<= 192.82 32 0.7 2.2 

HOC+3^2'^ 268.45 19 0.4 1.2 

^ S/N estimates do not include calibration error, which is 

estimated to be 10%. 

^ The unidentified feature at 290 GHz is described in ^T2\ 
° These two entries should be regarded as only upper limits. 

the starburst ring will also change. B would be expected 
to be for a beam-filling source and 2 for a point-source. 
We note that at these frequencies, the beam-scaled grey- 
body (first term) approximates to a powerlaw, where 

F^ocAflTi^^^ i^^+f^. (2) 

This means that our best-fit A value is really degenerate 
with n, and T, and B is degenerate with (3. Our best-fit 
model yields A = 0.0446 ±0.0002 and B = 0.421 ±0.027. 
Our goal with this model is to determine the best-fit con- 
tinuum in order to measure the line fiuxes accurately, not 
necessarily to properly model the physical conditions of 
the continuum. Therefore, our best-fit parameters may 
be degenerate with uncertainties in other physical pa- 
rameters as described above. For example, if we were 
to assume /3 < 2, the difference would be made up in 
the parameter B. A is much less than 1 likely because 
the dust continuum emission does not entirely fill the 
beam (geometrically), and the gas emitting the contin- 
uum likely has some filling factor less than 1. We would 
still expect B > because the changing size of the beam 
does encompass a varying amount of the starburst ring 
with frequency (i.e. we expect a frequency dependence 
for the relative fraction of coupling to the ring) . Finally, 
we note that the core/jet emission makes only a small 
contribution to the overall continuum, but it does affect 
the shape of the continuum and especially improves our 
fit at the higher end of our band; without it, we would 
likely overestimate the CS J = 6 — 5 flux by over 10%. 
In summary, this model adequately fits our continuum 
spectrum with reasonable physical parameters. 

The resulting measurements are given in Table [T] This 
fit had a reduced of 5.1 with 140 degrees of freedom. 
As can be seen in the bottom panel of Figure [U a large 
percentage of the total comes from the CO line, which 
is not well fit by a single-component Gaussian. We do not 
have the resolution to be able to better model this line, 
whether its assymetry is due to kinematic structure or 
blending of other lines. If we exclude the bins containing 
the CO line, the total is < 3. However, we do not 



use CO in our radiative transfer analysis, so its poor fit 
is not problematic for the rest of our results. 

Our fit includes one unidentified spectral feature at a 
rest frequency of approximately 290 GHz. This frequency 
may correspond to CH3CCH J = 17 ^ 16 or H2CO 
J = 4o4 — >■ 3o3, but we cannot identify other features 
of these molecules that might be in our band. Likely, 
this feature is a combination of emission features that 
we cannot resolve, and is therefore left as unidentified. 
Though other features may be present at the upper end of 
the spectrum, the channel-to-channel variation increases 
significantly above the CS J = 6— J-S transition due to im- 
perfect atmospheric subtraction, so we do not attempt 
to define spectral features above this frequency. SiO 
J = 7 — > 6 could possibly be identified in this range, but 
because we cannot identify the two lower-J transitions in 
our band, we do not fit this feature either. Below this fre- 
quency range, there arc a few other features of note that 
are not fitted. The first is a possible 3-sigma feature at 
approximately 284 GHz; given our inability to know the 
exact line center and no particularly strong transition ex- 
pected that at wavelength (though there are a few tran- 
sitions of methanol, CH3OH, around that frequency), we 
do not fit this line. Though the same difficulties apply 
to the unidentified line at 290 GHz, its significance was 
twice that of this line, and warranted inclusion in the fit. 
The apparent features at 200 and 260 GHz are not lines, 
but part of the continuum model, having been passed 
through our measured line profiles which included the 
instrumental sidelobes of the CO J = 2 — > 1 transition. 
Future use of Z-Spec will include modification of the line 
profiles to exclude these artifacts. 

2.3. Spectrum Results 

When comparing our fluxes to those previously pub- 
lished, care must be taken to account for different beam 
sizes because of beam dilution. Throughout this paper, 
we correct for beam dilution by dividing all Rayleigh- 
Jeans velocity integrated temperatures by 

n2 



where Oheam is the full-width at half maximum of the 
beam profile and Osource is the FWHM of the source. 
For all high density tracers, we assume a source size of 
4" (see discussion on the geometry of NGC 1068 in SJT]). 
Other s have been assuming smalle r source sizes such as 
1.5" ( Perez-Beaupuits et al.ll2009D . but we have found 
that a smaller source size does not significantly change 
our results for the physical conditions of the CND. 

For the rest of this paper, we will be focusing on high- 
density tracer emission from the CND, but first we find it 
appropriate to comment on the most prominent feature 
of our spectrum, ^^CO J = 2— > 1, mostly from the star- 
burst ring. A source size of 30" has been estimated from 
maps of CO emission (representing the diameter of the 
starburst ring, see After correction for source size 
and comparing to the one reported J = 2 — >■ 1 transition 
from iPerez-Beaupuits et al.l (120071) and three reported 
J = 2 — » 1 transitions from [Israeli ()2009l ) . Z-Spec's mea- 
sured flux is slightly higher than all four others. How- 
ever, after adding in a 15% calibration uncertainty to 
each measurement, they agree within the uncertainties. 
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This gives us an indication that our cahbration is consis- 
tent with others. 

The transitions useful to the study of the dense nuclear 
core of NGC 1068, the subject of the remainder of this 
paper, are those of HCO+, HCN, HNC, CN, and CS. 
The CS J = 4 — )■ 3 and CS J = 6 — ?► 5 measurements are 
new, and when combined with previous CS J = 3 — > 2, 
CS J = 5— >-4, and CS J = 7— >-6 measurements, provide a 
well-sampled ladder for radiative transfer modeling. We 
do not include CN in our models because we cannot ac- 
curately separate its hyperfine transitions; our reported 
CN intensity includes all 18 hyperfine lines included in 
CN J = 2^1. 

3. MODELING OF MOLECULAR GAS PHYSICAL 
CONDITIONS 

We sought to estimate the most probable physical con- 
ditions of NGC 1068's dense core by simultaneously mod- 
eling a number of molecular lines. Rather than just find 
the single most probable set of conditions, we examined 
a large parameter space and determined the most proba- 
ble set of conditions by marginalizing over all possibilites. 
To do so, we used B ayes' theorem to describe the likeli- 
hood of a particular set of model parameters g iven our 
observ ational measurements, as described in iWard et al.l 
(|2003i) . We have extended th is analysis to exami ne mul- 
tiple molecules at once, as in lNaylor et al.l (|2010D . Here, 
we first describe RADEX, which is used to create a grid of 
expected line fluxes for a range of physical parameters, 
then we detail the likelihood analysis and the specific 
considerations for NGC 1068. 

3.1. RADEX 

To investigate the physical parameters, such as tem- 
perature and density, that produced our measured line 
fluxes, we used RADEX, a non-LTE code freely avail- 
able from the Leiden At omic and Molecular Database 
(|van der Tak et all I2007D . We used the large velocity 
gradient (LVG) model to perform statistical equilibrium 
calculations using the escape probability in an expand- 
ing spherical cloud, but the results are insensitive to the 
precise form of the escape probability. When provided 
with the proper input data, RADEX calculates the level 
populations in the optically thin limit considering the 
background radiation field and then calculates the optical 
depths for the lines. The code continues calculating new 
level populations using new optical depth values until 
the two converge on a consistent solution. The program 
then outputs the resulting line intensities as background- 
subtracted Rayleigh- Jeans equivalent radiation temper- 
atures. 

The background radiation field can be as simple as 
a 2.73 K blackbody to represent the cosmic microwave 
background (CMB). However, the continuum emission 
from the CND should be used as the background if it 
dominates over the CMB; using the following procedure, 
we determined that this is not the case and that the CMB 
is an adequate background radiation field . We compared 
the core continuum measurements of Kri ps et al.l (|2006f ) 
to the CMB in our frequency range (88 to 363 GHz). 
They measure the continuum emission as a power-law 
with frequency due to optically thin synchrotron radi- 
ation from the central jet source, which is embedded 
within the CND but smaller than the full 4" diameter. 



TABLE 2 

RADEX Model Parameters and Ranges 



Parameter 


Range 


# of Points 


Tkin [K] 


100-7 _ 102-7 


45 


n{H2) [cm-3] 


102-5 - 10*-5 


45 


Ncs [cm-2] 


IQlO _ 1019 


40 


^HCO+/^CS 


10-1-5 - 101-5 


17 


^hcn/Xcs 


10-1-0 _ 102-0 


17 


Xhnc / ^cs 


10-1-5 - 101-5 


17 


Ay [km s-1] 


1.0 


fixed 


'-^background [K] 


2.73 


fixed 



Note. — All parameters are sampled evenly 
in log space. Ncs is given beca use it is defined 
as the primary species, see i|3.1l 



We use their measured emission sizes at 1 mm and 3 
mm to estimate a power-law relating source size to fre- 
quency to convert from flux density to specific intensity. 
However, the property which truly matters to the back- 
ground radiation field is the mean intensity, which is not 
equivalent to the specific intensity because the size of 
the continuum emission is smaller than the full extent of 
the 4" CND. Therefore, the mean intensity is dependent 
upon the angular extent of the central source, which is 
a function of the frequency-dependent size of the source 
and the distance in the CND from it. We use the same 
relationship between source size and frequency as men- 
tioned above and determine an average angular size of 
emission throughout the CND, which we use to calculate 
the mean intensity. We find that the mean intensity of 
the CMB radiation is comparable to that of our core -I- 
CMB model, and the difference in the predicted line in- 
tensities using either background model (CMB alone or 
core + CMB) is apparent in the RADEX results only at 
the lowest temperatures and densities; the overall likeli- 
hood results using grids with or without this background 
intensity added to the CMB are essentially indistinguish- 
able. However, we note that future studies using higher- 
J transition lines may still need to consider this back- 
ground, which dominates over the CMB above 500 GHz. 
In summary, the CMB -|- continuum emission is for the 
most part indistinguishable from the CMB alone, and so 
we use just the CMB for our background radiation field. 

The input parameters required for RADEX are black- 
body temperature of the background radiation (for our 
purposes, Tcmb — 2.73 K), kinetic temperature (Tkm), 
molecular hydrogen number density (assumed to be the 
only collision partner, n(7?2)), column density of the 
molecule (Nmoi), and the width of the molecular lines. 
Our calculations are all computed per unit linewidth be- 
cause this is the physically relevant quantity that deter- 
mines the optical depth scale; we later scale this value 
by the linewidth in our likelihood analysis. Furthermore, 
RADEX assumes that the emission region fills the entire 
beam, but our likelihood analysis also introduces a beam 
filling factor <I>^ to account for dumpiness in the gas. 

We used RADEX to create a grid of line intensities for 
a range of Tkin, n{H2), and Nmoi- RADEX only treats 
one molecular species at a time, therefore in order to con- 
duct a multi-species analysis, we first created a RADEX 
grid for one species, the primary species (CS). We next 
created grids for all other (secondary) species, in which 
we also explore their relative abundances to the primary 
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species. In calls to RADEX, the column density was 
the product of the primary species' column density and 
the secondary species' relative abundance. Our analysis 
used a relatively coarsely sampled but wide ranging grid, 
detailed in Table El Though each RADEX grid is cre- 
ated separately, the presence of multiple species at once 
is simultaneously considered in the likelihood analysis. 

3.2. Bayesian Likelihood Analysis 

We next compare the calculated line intensities to the 
measurements, detailed in Table |3l We add calibra- 
tion error to the measurements in quadrature with the 
line uncertainties, assuming independent Gaussian un- 
certainties; this assumption is discussed in t j4.2l 

Given a set of measurements x and model parameters 
P = {Ncs,n{H2),Tkin,^A,Xjnoi/Xcs), the Bayesian 
likelihood of the model parameters given the measure- 
ments is 



P{p\x) = 



P{p)P{x\p) 
P{x) 



(4) 



where P{p) is the prior probability of the model parame- 
ters (see §3.3|) . P{x) normalizes, and P{x\p) is the prob- 
ability of obtaining the observed data set given that the 
source follows the model described by p, which is the 
product of Gaussian distributions in each observation, 



P{x\p) = TT — ^ exp 



2a t 



(5) 



Above, Ui is the standard deviation of the observational 
measurement for transition i and /^(p) is the RADEX- 
predicted line intensity for that transition and model. 
The product is carried out over all of the molecular 
species considered. 

3.3. Prior Probabilities 

The prior probability allows us to include previously 
known information about NGC 1068 in our likelihood 
analysis. We created a "binary" prior in which all phys- 
ical situations were assigned P{p) — 1, and all geomet- 
rically unphysical situations were assigned P{p) = 0. 
There were 4 criteria that each point in our model grid 
needed to satisfy in order to be deemed physically plaus- 
able. 

Tau — At high optical depths, RADEX is unreliable be- 
cause the cloud excitation temperature can become too 
dependent on optical depth with large column densities. 
Therefore, we only include results with r < 100 as rec- 
ommended by the RADEX documentation. Most of the 
grid points excluded by this prior are at very high column 
density; in fact, the likelihood results are nearly identical 
with or without this prior condition. The average opti- 
cal depths for the columns that we find for our likelihood 
solutions are generally less than 10. 

Total mass — The mass in the emission region {Mregion) 
should not exceed the dynamical mass of the galaxy, so 
we require that 



1 



< M, 



dyn ■ 



(6) 



We are assuming the emission region is 312 pc wide (see 
21, so Aregion is the Corresponding area. Nmoi^A is the 
beam-averaged column density, m//^ is the mass of the 
hydrogen molecule, the factor of 1.5 accounts for helium 
and other heavy elements, and Xmoi is the abundance of 
t he molecule r e lative to hydrogen, rimoi /nH2 ■ 

lUsero et"ari (|2004[ ) conducted previous LVG simula- 
tions to determine the relative chemical abundances of 
many of the molecules discussed here; they found Xcs 
to be 2.0 X 10"*^ in the East knot and 1.6 x 10"** in 
the West knot. These values were derived assuming 
Xco = 8 X 10~^. Our beam does not distinguish be- 
tween the two regions, so as a conservative limit we use 
the higher value. For any grid point of a given column 
density and filling factor, using a higher value of X^oi 
results in a lower estimate of the mass, meaning that 
fewer grid points overall will exceed the dynamical mass 
cutoff. 

The dynamical mass can be estimated as 



M, 



dyn 



G ' 



(7) 



where cr^ is the projected (ID) stellar velocity dispersion 
over an aperture of radius R, and the 2 accounts for 2D 
motion. lOliva et all (|1995l ) measured a velocity disper- 
sion of 161 ± 20 kms~^ in a 4.4" beam, which closely 
approximates the size of the region which we are study- 
ing. This yields a dynamical mass of 2 xlO^ Mq, which 
we use as an upper limit for the mass constraint. 

Velocity Gradient — We can estimate the observed veloc- 
ity gradient of individual molecular clouds, and require 
that it be at least the minimum gradient required for 
virialization. Using the expressions for these two quan- 
tities. 



cs 



Ncs 



and 



= [ -iraGnmH^nH^ 



1/2 



(8) 



(9) 



where a = 1, = 1.5, and mjj^ is the mass of the hy- 
drogen molecule, we can place the requirement that 



(dv/dr)obs 
(dv/dr)vir 



> 1. 



(10) 



Column length — Finally, the length of the column {Lcoi) 
of the primary species should not exceed the length of the 
entire molecular region. We assume that the length in 
the plane of the sky is the same as that orthogonal to the 
plane of the sky. For the 4" region we round the length 
of the region to 350 pc for the purposes of creating an 
upper limit. 



N„ 



< L 



(11) 



We find that the length prior has the greatest effect on 
the likelihood results by excluding the lowest molecu- 
lar hydrogen density and highest column density solu- 
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tions. Though it is the combination of these two param- 
eters that matters in the length prior, the final marginal- 
ized likehhood resuhs generally exclude solutions below 
^ IQ3.5 QT[^-3 (iensity and above ~ 10^^-^ cm~^ in CS 
column density that would otherwise be allowed without 
this prior. 

In summary, the prior probability is defined as: 



'1 ; T < 100, M, 



region 



< M, 



dyn-) 



; otherwise 



(12) 



The parameters that are assumed for the likelihood anal- 
ysis are summarized in Tabled 

3.4. Contamination of CND intensity by Starburst Ring 

iPerez-Beaupuits et al.l ()2009f ) noted the ~ 30" diam- 
eter starburst ring could be contributing high-density 
tracer flux measureable by a large enough beam; be- 
cause we are only interested in modeling the CND of 
NGC 1068, this extra fiux would need to be removed 
from our calculations. Though other single-dish mea- 
surements have sufficiently small beam sizes to avoid this 
problem at higher-J transitions, Z-Spec's beam is large 
enough (26-37" FWHM) that we do need to take this into 
consideration. They constructed a first-order estimate 
(their equations 1-4) of what they call the starburst con- 
tribution factor {/sb) by comparing two measurements 
of the same line with different sized telescope beams. 
The fiux emergent from only the CND would be the to- 
tal flux measured times (l-fss)- This estimate, however, 
yields high starburst contribution factors (around 50%, 
with large uncertainties), which disagrees with currently 
available interferometric observations (such as those ref- 
erenced in 21) • 

We attempted to calculate fsB for the Z-Spec mea- 
sured lines by comparing them to other measurements in 
the scientific literature. In this paragraph, in parenthe- 
ses after each citation, we quote the reported velocity- 
integrated line intensities in Kkms~^ and conversion to 
Jy kms^"'^ based on the references' reported beam sizes 
(the total integrated fiux density in Jy kms~^ is the 
appropriate quantity to compare between telescopes j3- 
The Z-Spec measurement to which we are comparing 
is reported in Table [TJ In our calculations, we add 
20% calibration error into the line measurements (this 
is slightly more than used elsewhere in the paper be- 
cause of the additional uncertainty of the assumptions 
that went into the estimate itself, such as Gaussianity 
of the beam). For the HNC J = 3 — >■ 2 line, we com- 
pare to iPerez-Beaupuits et all (|2007D (3.5 Kkms'S 69 
Jy kms~^) and find that for our beam, fsB = 0.54±0.37 
(0.30 ± 0.31 for th eir JCMT beam). In HCO+ J = 3->2 
we can compare to lKrips et al.l ()2008[ ) (7.6 Kkms^^, 40 
Jy kms~^), who measured this transition using a 9.5" 
beam. At such a small beam size, the first order esti- 
mate essentially forces our measurement to match theirs, 
with fsB = 0.85 ± 0.05. HON J = 3 ^ 2 was measured 

''' Note that this con versi on is different than the beam-dilution 
correction described in ^2.31 which assumes a source size (and pro- 
duces the Beam Corrected Intensity column of Table [3}. In this 
section we are simply converting to flux density units. 



by IKrips et all (l200l) (19.0 Kkms-\ 99 Jy kms-^), 
IPerez-Beaupuits et all (pOOTi) ( 22.7 Kkms-\ 425 Jy 
kms-i), and iBussmann et al.l 12008") (6.01 Kkms-\ 
313 Jy kms~^) but we find inconsistencies when com- 
paring to these three other observations. The last two 
fluxes in Jy kms~^ are comparable to ours despite dif- 
ferent beam sizes (implying that we are measuring the 
same emission, which must be concentrated in the cen- 
ter if all beams are measuring it); this yields a negative, 
non-sensical estimate of the flux from the starburst re- 
gion. We can co mpare to the smaller (9.5") beam of 
iKrips et al.l (|2008( ) again yielding our fsB = 0.76 ± 0.07. 
CS yields 3 measured transitions in our band, only 1 of 
which ( J = 5 — >■ 4) we are able to compare to other mea- 
surements from iMiltlneFan (HOOl) (3.3 Kkms^i, 16.2 
Jy kms~^) with a beamsize of 10". Our fsB for this line 
is 0.72 ±0.1. 

These comparisons yielded suspiciously high contribu- 
tion factors from the starburst ring. The fact that our 
measured fs b values for the J = 3 — >■ 2 transitions with 
me dium-sized beams (29" ) are h igher than those derived 
by iPerez-Beaupuits et al.l ()2009() for J = 1 — > and much 
larg er beams (44 - 55") adds to this suspicion. Addition- 
ally, lUsero et al.l (j2004[ ) also made estimates of the con- 
tamination from the starburst ring by comparing central 
pointing values to those offset from it and found little 
contamination for SiO J = 3— ?>2, no more than 25% con- 
tribution for SiO J = 2 1, and no more than 30% for 
II^'^CO+ and IICO+ J = l^-0. Furthermore, in the case 
of HCN J = 3 — 2, we are already consistent with two 
other measurements cited above without any correction 
factor. A likely explanation is that the values to which 
are we comparing are too low, potentially due to errors 
in calibration or baseline continuum measurements. An- 
other possibility is that the azimuthal structure of the 
starburst ring (or the beams) is significantly affecting 
the calculation of /sb, which assumes such symmetry. 

If we were to correct our measurements using the star- 
burst correction factor described above, we find that the 
resulting spectral energy distributions of the CND for 
these high density tracers cannot reasonably be described 
by a single-tcmpcraturc component of gas given the re- 
ported J 4 ^ 3 intensities of HCO+, HCN, and HNC 
(jPerez-Beaupuits et al.ll2009l) . which rise above our cor- 
rected J = 3 — >■ 2 intensities when corrected for source size 
(as described in ij2.3p . There are two possibilities: 1) the 
J = 4— ^3 lines are primarily excited in a separate compo- 
nent of gas within the CND that is warmer than the gas 
primarily traced by the lower J lines, or 2) the J = 3— >2 
lines to which we are comparing in order to derive fsB 
are too low, perhaps due to pointing or calibration errors. 

Because we currently have an incomplete picture of the 
extent of the contribution to these line strengths from the 
starburst ring, we present results for 3 scenarios. The 
first is our preferred model, whose results will be pre- 
sented in full. We only make one correction for starburst 
emission in this model. The majority of the lines used 
in our analysis have critical densities higher than SiO 
J = 3 2 {nc„t = 9.6 X 10^ cm-^ at 60 K), with the 
only exception being HCO+ J = l— >0. Based on Usero's 
measurements (which found that SiO J = 3 2 had no 
contribution from the starburst ring), we assume in our 
preferred model that all of our molecular transitions are 
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TABLE 3 

Observed Line Intensities Used in Likelihood Analysis 



















Beam 


Beam Corrected 




Species 


Transition 


I^rest 






Intensity 


FWHM 


Intensity'' 


Ref 








[GHz] 


[K] 




[Kkms-i] 


["] 


[Kkms-i] 




CS 


J 


= 3 -5> 2 


146.97 


14.11 


1.3 


xlO^ 


9.5 ± 1.4 


16.0 


162 ± 34 


1 






4^3 


195.95 


23.51 


3.1 


xlO'' 


1.56 ± 0.43 


37.1 


136 ± 40 


2 






5^-4 


244.94 


35.27 


6.9 


xlO® 


1.23 ± 0.23 


30.8 


74 ± 16 


2 






6 -5> 5 


293.91 


49.37 


1.2 


xlO'^ 


1.17 ± 0.39 


26.3 


52 ± 18 


2 






7^6 


342.88 


65.83 


2.0 


xlO'^ 


1.4 ± 0.5 


14.0 


18.6 ± 7.2 


3 


HCO+ 


J 


= 1^0 


89.19 


4.28 


2.2 


xlO^ 


14.6=^ ± 0.2 


29.5 


809 ± 122 


4 






3^2 


267.56 


25.68 


3.8 


xlO^ 


5.58 ± 0.35 


29 


299 ± 35 


2 






4^-3 


356.73 


42.80 


9.1 


xlO'' 


3.8 ± 0.5 


14 


50 ± 10 


5 


HCN 


J 


= 1 -> 


88.63 


4.25 


3.2 


xlO" 


10 ± 1.7 


44 


1220 ± 177 


6 






3^2 


265.89 


25.52 


5.2 


xlO'' 


8.5 ±0.26 


29 


455 ± 48 


2 






4 -> 3 


354.51 


42.53 


1.2 


xlO* 


13.9 ± 1.6 


14 


184 ± 35 


5 


HNC 


J 


= 1 -> 


90.66 


4.35 


3.6 


xlO" 


11.4 ±0.7 


25 


457 ± 74 


7 






3 -5> 2 


271.98 


26.11 


5.8 


xlO'^ 


2.15 ±0.21 


28 


108 ± 15 


2 






4 -5> 3 


362.63 


43.51 


1.4 


xlO* 


2.7 ±0.3 


14 


35.8 ± 6.7 


5 



References. — (1) [Mauersberger et al. (IgS^), (2) Z- Spec (tliis work) . f 3) IBavet et aLI j200g ) . 

(4) iK rips et a l.i 12003'), (5) Pcrez-Beaupuits et al] 120091 '). (6) 'Pcrez-Bea upuits et al.l 12007l l. T?) 
IHuettemeister et al.l ljl995l '). 

Calculat ed at 60 K for HON and HNC, 70 K for CS and HCO+from the Leiden Atomic and Molecular 
Database fScho ier et al.|[2005l ') ■ where n^TU = Aji;/7„;. 

^ See i|2.3l for details. Also includes 10% calibration errors for Z-S pcc and 15% for others. 

Later multiplied by 0.7 ± 0.1 for the Preferred model, see i|3.4l for details. Other lines are corrected 
for other models by factors described in the aforementioned section of the text. 



TABLE 4 
Likelihood Parameters Used 



Parameter Value Units 



Line width'' 


240 


[km s-i] 


Abundance {Xcs/^H^)^ 


2.0 xlO"* 




Angular Size Scale'^ 


78 


[pc/"] 


Source size'^ 


4.0 


["] 


Length Limit*^ 


350 


[pc] 


Dynamical Mass Limit'' 


2 xlO^ 


[Mq] 



" Used for scaling the line intensities, see ^2.21 
^ Used for prior probabilities, see i|3.3l Relative 
abundance is for the primary species, CS, relative 
to Ha. 

Used for prior probabilities (source size also used 
for scaling to correct for beam dilution) , see [[l] 

only excited in the CND, except for HCO+ J = 1 -J' 
which we decrease by 30%. This imphes that the gas in 
the CND is denser than in the starburst ring. "We beheve 
this is the best model at the current time, though future 
interferometric studies may present new information that 
modifies this scenario. 

The next two scenarios are presented in the Figures [H 
[3l m [5l and El but not in any tables, though we discuss 
the differences in SHI The scenario labeled "1-0 /sb" 
uses the fsB values of Pcre z-Beaupuits et al.l ([2009) for 
the HCN, HCO+, and HNC J = 1 ^ values (/ss = 
0.56, 0.45, and 0.45, respectively). The final scenario la- 
beled "/sB, Cool" uses the same fsB values as in "1-0 
/sB," but also uses the fsB for Z-Spcc measurements as 
described earlier in this section (for HCN, HCO+, and 
HNC J = 3^2 and CS J = 5^4, Jsb = 0.76, 0.85, 0.54, 
and 0.72, respectively) Because this reduces the J = 3— >-2 
lines, we exclude the J = 4— >-3 lines, assuming they may 
be tracing a separate, warmer component. We also in- 
clude one extra line in this model, HCN J = 2 — > 1 from 
iKrips et"al., (,2008.) . which had been excluded from the 



other scenarios because of its low value. The spectral 
energy distributions in 2] should clarify the cases pre- 
sented. To summarize, the preferred model reduces only 
the HCO+ J = 1 ^ line by 30% due to starburst con- 
tribution, "1-0 /sb" corresponds to reducing emission 
of the ground transitions of HCN, HNC, and HCO+ by 
about half (assuming the other half is actually emission 
from the starburst ring), and "/sb. Cool" additionally 
corrects for a large starburst contribution in the middle- J 
states and excludes higher-J states which may be tracing 
a warmer component. 

4. MODELING RESULTS AND DISCUSSION 

The raw result of our analysis is a matrix of likelihoods 
each characterized by a certain value of T^im '"■(-^2), 

Ncs, ^Ai Xhco+ I ^cs-, Xhcn / Xcsi and Xhnc/Xcs- 
To obtain a distribution for an individual parameter, 
we integrate the likelihood matrix over all other dimen- 
sions. The integrated distributions for the first four pa- 
rameters are shown in Figure [2J We can use these pa- 
rameters to create distributions for 3 other supplemen- 
tary parameters (bottom of Figure [3]) : thermal gas pres- 
sure, P — X Tkin, beam-averaged column density, 
< Ncs >= X Ncs, and velocity gradient of individ- 
ual clouds, Equation[51 From the beam-averaged column 
density we calculate the total molecular mass in the beam 
with Equation [BJ which is shown in the upper x-axis in 
the bottom middle panel of Figure [H 

The column densities for secondary molecules are 
shown in Figure |4l The results for the distributions, 
marginalized over all other parameters, are in Table [S] 
(for the "Preferred" model only). The resulting spectral 
energy distributions and their optical depths are shown 
in Figures [5] and [SI Uncertainties cited in text are la 
unless otherwise stated. 

4.1. Physical Conditions 
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TABLE 5 

Multi-Species Radiative Transfer Results of Preferred Model 



Parameter 


Median 


la Range 




Max 


4D Max 


Tkin [K] 


117 


44 


- 231 




159 


159 


n(H2) [cm-3] 


3.4 X 10* 


1.7 X 10* 


- 8.7 X 


10* 


3.5 X 10* 


2.6 X 10* 


Ncs [cm-2] 


1.9 X 10^5 


1.0 X 10^5 


- 3.7 X 


10^5 


2.4 X 1015 


2.4 X 1015 


*A 


0.48 


0.30 


- 0.72 




0.49 


0.62 


P [K cm~^] 


4.0 X 10" 


2.3 X lO'' 


- 6.9 X 


10'= 


4.7 X 10" 


4.7 X 10" 


< Ncs > [cm-2] 


9.3 X IQi* 


4.2 X lOi* 


- 1.6 X 


1015 


1.1 X 1015 


1.1 X 1015 


dv/dr [km s"-*- pc"-*-] 


167 


51 - 


- 581 




212 


212 


Mass in Beam [Mq] 


8.5 X 10^ 


3.8 X 10^ 


- 1.4 X 


10« 


1.0 X 10* 


1.0 X 10** 


^hcn/^cs 


4.8 


3.1 


- 7.6 




4.9 


4.9 


-/Vhcn [cm-2] 


7.9 X IQis 


2.8 X 10^5 


- 2.3 X 


10i« 


2.0 X 10" 


2.0 X lOis 


^HCO+ / ^CS 


0.30 


0.23 


- 0.47 




0.42 


0.42 


^^HCO+ [cm"^] 


7.2 X 10^"' 


2.6 X 10^* 


- 1.4 X 


1015 


8.9 X lOi* 


8.9 X IQi* 


JfHNc/-''^CS 


1.07 


0.79 


- 1.66 




1.55 


1.55 


TVhnc [cm-2] 


2.0 X 10^5 


6.9 X 10^* 


- 5.5 X 


10l5 


2.4 X 10i5 


2.4 X 10i5 



Note. — ID max refers to the maximum value of the integrated parameter distri- 
bution. 4D max refers to the value of that parameter at the best fit solution. 



TABLE 6 
Relative Abundances 



Abundance Ratio 


Median Value 


^hcn/^hco+ 


15.8 ± 10 


^HNC/-''^HCN 


0.22 ±0.15 


-'<'CS/''<^HCN 


0.21 ±0.1 


-''^hcn/^co 


0.0012 ±0.0005 




2.5 X 10"* 



Note. — All values are derived 
from Table \5\ CS and CO use as- 
sumed abundances to H2 of 2 X 10~* 
and 8 X 10~5^ respectively. Uncertain- 
ties are representative. 

Warm gas (44 — 231 K) is implied in the nuclear disk. 
Because the temperature distribution is so broad, it is 
possible that a single-temperature model is not adequate 
for the CND of NGC 1068. It is important to note that 
the temperature and the density are degenerate; a high- 
temperature, low-density model may reproduce the data 
equally as well as a low-temperature, high-density model. 
Therefore, it is not surprising that both our temperature 
and density distributions are somewhat broad. How- 
ever, for the preferred model, low temperatures (< 50 
K) are disfavored, which corresponds to the upper limit 
on the density likelihood distribution. The density is 
constrained to lO**'^ — lO*^'^ and the product of the tem- 
perature and density (pressure. Figure [3|) is better con- 
strained with Log(P [K cm~3]) ^ q qq ^ o.24). We note 
that decreasing the HCO"'' J = 1 flux can bias the 
results towards higher pressures (as can be seen in the 
higher pressures of the other models, bottom left panel 
of Figure|3|); however, as will be discussed in M.2\ HCO+ 
does not dominate the results. 

Previous models have assumed temperatures of around 
50 K or 80 K, but we have shown that the gas (or at least 
a significant component of it) likely traces a higher tem- 
perature (^ 159 K) region. This distinction matters be- 
cause assuming a lower temperature necessarily, and per- 
haps incorrectl y, increases the den sity that is required to 
find solutions. iKrips et"aLl ()2008[ ) found 2 best fit mod- 
els for NGC 1068, one low-temperature/high-density and 
one high(unconstrained)-temperature/low density. By 




100 200 300 400 
T.„ [K] 
Log,„ N[H,^ [cm"'] 



3 4 5 6 7 
Log,o n[H3] [cm"^] 



^ 0.2 




15 16 17 
Log,, N[CS] [cm"'] 



Fig. 2. — Multiple-species probability distributions for the four 
primary parameters, T^in, n^^, Ncs> a-nd All are normalized 
such that the peak maximum likelihood value is equal to 1. The 
hydrogen column density axis uses our assumed CS abundanc e of 
2 X 10~* relative to H2. The different models as described in i|3.4l 
are displayed with different linestyles: solid ("Preferred"), dashed 
("1-0 /sfl"), and dash-dotted ("/gs, Cool"). 



examining the relative likelihoods over a large parameter 
space, we have demonstrated that there is broader sup- 
port for their second solution, which matches well with 
our measurements. At the least, it is clear that a high- 
temperature component is needed to adequately explain 
the measured spectral line energy distributions. 

The column densities and relative abundances of the 
molecules are well constrained, as can be seen in Figure 
|3]and Table [5] Moreover, the column densities are fairly 
consistent among the different models, though CS and 
HNC "switch" places relative to one another in bottom 
two models/panels in Figure |4j however, their distribu- 
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Fig. 3. — Multiple-species likelihood results: 2D distributions and supplementary parameters. The top row contains 2D probability 
distribution contours for the parameters from Figure [2] each contour level represents 0.8, 0.6, 0.4, and 0.2 of the maximum of that 
distribution. Diagonal lines on the contour plots correspond to values of the bottom para mete r (P, N(CS), and dv/dr from left to right), 
which are derived from the primary parameters. The different models as described in il3.4l are displayed with different colors (top) or 
linestyles (bottom): black/solid ("Preferred"), blue/dashed ("1-0 /sb"), and red/dash-dotted ("/ss. Cool"). 



tions still overlap in all models. Another strongly con- 
strained parameter from these models is the total mass 
in the beam, bottom middle panel of Figure [3l at 10* 

Mq. 

We have referred all of our measurements to a 4" source 
size and find, with that referred size, an area filling factor 
of 0.48tQ 28 (for the preferred model), which corresponds 
to a volume filling factor O.SSIq^^ (assuming spherical 

symmetry and therefore $y — ^'^^)- This implies some 
dumpiness in the gas. Our infer red mass in t his 4" area 
is 1.0 xlO* Mq, very similar to [Israeli (|2009[ )'s estimate 
of 1.2 xlO* Mgas inferred from CO. 

The derived velocity gradient is much greater than 
what may be estimated by simply dividing the total ob- 
served velocity linewidth by the length of the CND (~ 
1 kms^^pc^^). This indicates that the CND is likely 
composed of clumpy, turbulent structures. Such high ve- 
locity gradients have been found in studies of the gas 
in the NGC 253 starburst nuc leus ()Bradford et al.ll2003l: 
IHailev-Dunsheath et al.l I2008D and can arise naturally 
due to the geometry and k inematics in compact circum- 
nuclear environments (e.g. [Bradford et al.l (|2005[ )). 

In our preferred model, we have only corrected HCO"*" 
J = 1 — > for contributions from the starforming ring (see 
ij3.4l for justification, and ij4.2l for discussion of the effects 
of this choice). If we use th e Jsb correction factors of 
iPercz-BeaupuitseLaLj (|2009[ ) ("1-0 /ss", dashed line in 
Figures [2] and 13]), as well as those for the J = 3— >-2 lines 



that we can derive by comparing Z-Spec's measurements 
to other data ("/ss. Cool" , dashed-dotted line in Figures 
|2]and|3]), the likelihood results (especially for tempera- 
ture and density) are dramatically different. The temper- 
ature peaks at 5 K, but is unconstrained (constant rela- 
tive likelihood) at higher temperatures. As a result, the 
density distribution is significantly broader and uncon- 
strained even above 10^ cm~^. The resulting SEDs (Fig- 
ure [5]) , after likely overcorrecting for contribution from 
the starburst ring, do not seem to be describing physi- 
cally reasonable parameters. Even if we only use their 
derived J = 1 — correction factors and leave all other 
measurements uncorrected, the temperature still peaks 
at a very low value (12 K) and is unconstrained along 
with the density. However, as mentioned, the column 
densities are still fairly well constrained, and given this 
constraint the models show a trade-off between $^ (de- 
creasing with each model) and density (increasing with 
each model) . We believe the SEDs that we are modeling, 
where only HCO+ J = 1 — > needs a correction for contri- 
bution from the ring, are the most reasonable considering 
interfcromctric maps described in fJT] and also produce 
more physically reasonable results than if we apply extra 
corrections. However, we recognize that there may still 
be some level of contamination from the starburst ring 
that we cannot properly calculate, and therefore the re- 
sults may represent a combination of the CND and ring, 
but should be dominated by the CND. The effects of 
the assumed fs b values are further discussed in the next 
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Column Density [cm '] 

Fig. 4. — Multiple-species column density like lihoods. Each panel 
corresponds to one model as described in >13.4I 

section. 

4.2. Systematic Effects 

It is important to note a few systematic effects that 
exist within our likchhood analysis. We have assumed 
that the measurements all have independent Gaussian 
uncertainties. Such a treatment would only be com- 
pletely valid in the case that the calibrations for each 
line measurement were independent; however, the Z- 
Spec line fluxes have correlated calibration uncertain- 
ties. Additionally, multiple lines come from the IRAM 
30-meter telescope and the James Clerk Maxwell Tele- 
scope (JCMT), and the lines from each individual tele- 
scope have correlated calibration uncertainties. Because 
it would be difficult to assign relative calibration errors 
to the various line measurements (they are not reported 
uniformly), we do not attempt to treat the calibration 
errors as a separate systematic error. Rather, we simply 
assess the impact of the inclusion of the assumed rel- 
ative calibration errors in quadrature with the random 
errors. Doubling the calibration errors results in likeli- 
hood distributions that widen (e.g. by a factor of 1.4 
for pressure at the 1-sigma range), but the median val- 
ues change insignificantly. This indicates that line fluxes 
determine the molecular level populations and the sizes 
of the assumed relative calibration errors contribute to 
the likelihood widths. (Of note, larger calibration errors 
reduce the jaggedness of the temperature likelihood dis- 
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Fig. 5. — Line spectral energy distributions. Squares represent 
data points (see Table [3] for references used, filled squares indicate 
Z-Spec data points), X's represent the best fit model results. All 
measurements are scaled for beam dilution, RADEX results are 
scaled by linewidth. 15% calibration is added in quadrature for all 
measurements, except for Z-Spec, fo r wh ich 10% is added. Each 
column is one model as described in i|3.4l 

tribution in Figure [2] This implies that the grid point 
spacing is slightly sparse relative to the measurement un- 
certainties, but since the mean and overall shape do not 
change, there are no major consequences for our conclu- 
sions.) 

We use the linewidth to scale the intensities by first di- 
viding each velocity-integrated intensity by the assumed 
linewidth to compare directly to the RADEX predictions 
of Raleigh- Jeans temperature for a given column density 
per unit linewidth. Then, in order to plot the likelihood 
for the total column density, we multiply column den- 
sity per unit linewidth by the linewidth. We assumed 
the same linewidth (240 kms^^) for each molecul ar line 
transi tion based on the CO J = 2 — > 1 line of llsraell 
(|2009D . though there are slight variations by molecular 
species. For exa mple. Tables 3 of Pcrcz-Beau puits et al.l 
(2007) and 2 of iPerez-Beaupuits et al.i (|2009i) illustrate 
the linewidth for some of the molecular transitions used 
here; many are broken in multiple components, and some 
single-component measurements have linewidths as high 
as 275 kms^i (HCN J = 3^ 2) and as low as 166 kms'^ 
(HCO+ J = 4 — > 3). Changing the assumed linewidth 
modifies the likelihood results for the column density. 
Using a linewidth of 166 instead of 240 kms~^ reduces 
median CS column density by 25%, and using a linewidth 
of 270 instead of 240 kms~^ increased the median CS 
column density by less than 1% (the distributions for 
the area filling factor shift in the opposite manner from 
column density in similar proportions). Additionally, the 
supplemental parameter of velocity gradient would be af- 
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one model as described in N3.4I Optical depths are for the best-fit 
model results, as indicated in Figure [5] 



fected by using a different linewidth (in proportion with 
that hncwidth, i.e. dv/dr would be 30% smaller if the 
linewidth were 30% smaller). Thus, given this range of 
possible values, the systematic errors due to linewidth 
assumptions are less than 30%. 

The CS abundance was assumed; this does not affect 
the relative abundances of the secondary to the primary 
species, but would affect the conversion to molecular hy- 
drogen. Both the total mass in the beam and the velocity 
gradient depend (in inverse ways) on X^g. The fa c t that 
our estimate for the gas mass agrees with Hsraeil pOOl ) 
suggests that the value we have adopted for Xcs is about 
right. Finally, Xc^ is used in the prior probabilities to 
determine if individual grid points violate the maximum 
dynamical mass, length, or velocity gradient criteria. If 
Xc5 is too high/low then more/fewer grid points are 
allowed. The main effect of the prior probability is to 
cut the lowest density and highest column density points 
from consideration in our analysis. 

This analysis necessarily assumes that all the molecular 
lines trace the same gas. Detailed line profiles indicate 
that not all li nes necessarily trace th e sam e kinematic 
structures (i.e. iPerez-Beaupuits et al.l (|2007[ )): however, 
this analysis is an attempt to examine the bulk prop- 
erties of the CND. We investigated this claim by also 



examining the likelihood distributions created by mod- 
eling individual molecules one at a time. These single- 
species models are less constrained than the multi-species 
model (as one would expect by using fewer molecular 
lines), but show some mutually consistent results. The 
temperature and density distributions of HNC are most 
similar to the multi-species models. HCN and CS show 
similar profiles, and though their best-fit solution re- 
quires temperatures of ^ 160 and 80 K, respectively, the 
marginalized distribution also shows a pronounced peak 
at lower (^ 10 K) temperatures. Because of these allowed 
lower-temperature solutions, the density distribution for 
these two molecules peaks slightly higher than the multi- 
species model and is less constrained at the high-density 
end, though their uncertainties are quite large; the la 
uncertainties for density of HCN and CS span 1.5 and 
2.4 orders of magnitude, respectively. HCO+is the 
only molecule that, when modeled alone (including the 
30% correction to J = 1 — > 0), only demonstrates low- 
temperature/high-dcnsity solutions, though the density 
is also very unconstrained, as the la range spans 2.4 or- 
ders of magnitude. When we model only HNC, HCN, 
and CS (excluding HCO"*"), the results do not change 
significantly compared to when we include HCO^, indi- 
cating that is not dominating the multi-species likelihood 
results, but may be tracing a separate component of gas. 

Though we acknowledge that our treatment of the 
CND as one chemical/physical component is a simpli- 
fying assumption that is likely not strictly correct, our 
results still demonstrate fairly consistent results across 
multiple molecules. We present these results with the 
understanding that within this large region, there are 
likely multiple chemical/physical components, but with 
single-dish measurements we are only capable of model- 
ing and commenting on the bulk properties of the region. 
These physical parameters describe the sum of the emis- 
sion we see, but not the individual components. ALMA 
will be key to resolving individual components within the 
CND. 

Finally, we note that there are admitted uncertainties 
in the starburst contribution factors (fsB) for our pre- 
ferred model, which we have only corrected for HCO^ 
J = 1 — > 0, using fsB ^ 30%. This correction serves to 
slightly increase the density while decreasing the tem- 
perature and filling factor, though the column density is 
not significantly changed. The pressure distribution is 
only slightly increased (our model peaks at 4.7 x 10^ K 
cm^^, though it would peak at 3.1 x 10^ K cm^^ with- 
out the correction). The mass distribution is slightly 
changed as well, as without the correction it would peak 
at 2 X 10^ Mq, a factor of 2 higher than our model. In 
total, the HCO"*" J = l— )-0 correction changes the distri- 
bution only a small amount relative to the uncertainties. 
For the other models, we first note that the current val- 
ues of fsB used are subject to systematic errors in their 
derivation, because the model is sensitive to both the 
beam profiles, which may not be purely axisymmetric 
Gaussians, and also to the relative calibration between 
the telescopes (which is notoriously difficult). For any 
individual molecule, as we increase fsB (reducing the 
J = 1^0 intensity for HCN, HNC, and HCO+), the best- 
fit solution's density increases and is less constrained, and 
the relative likelihood of higher temperature decreases. 



DENSE MOLECULAR GAS IN THE CND OF NGC 1068 



13 



In general, the filling factor also decreases. Higher fsB 
scenarios require solutions which peak at the J = 2 — >■ 1 
line instead of J = 1 — )■ (see the middle column of Figure 
[5]) , which we cannot constrain given no J = 2 — >■ 1 mea- 
surements. Our slightly lower-density solutions may be a 
result of uncertainties in fsB, though we have attempted 
to model the CND given all the currently available infor- 
mation that we have. In the future, when fsB is bet- 
ter constrained through interferometric measurements, 
our models may require modification that will result in 
slightly higher-density solutions. 

5. DIAGNOSTICS OF THE AGN CENTRAL ENGINE 

The results of our radiative transfer analysis provide 
diagnostic information about the energetics of the CND. 
We first present a discussion of our derived relative 
molecular abundances and how they fit into the con- 
text of previous studies of the CND, especially its XDR 
nature. We next present a separate analysis of line ra- 
tios specific to the Z-Spec band, which also supports the 
XDR scenario. Finally, we speculate on the lifetime of 
the CND as determined by the black hole accretion rate 
and inferred molecular gas reservoir. 

5.1. Molecular Abundances and Possible Energy Sources 
for the Molecular Gas Excitation 

Some relative abundances of interest are presented in 
Table [51 We particularly note the high relative abun- 
dance of HCN to HCO"*". The median values from Ta- 
bleia yield [HCN]/ [HC0+1 of « 16. This is higher than 
lUsero et al.T(|2004l) (with ratios of 1.3 and 0.6 for the East 
and West k n ots, r espectively) but more consistent with 
IKrips et all (|2008l) . who required [HCN]/[HCO+] of at 
least 10 in order to find a solution that matched mea- 
sured values. 

There are three ways to produce such a high abun- 
dance. HCN could be enhanced in a PDR by far-UV 
radiation from massive star formation. However, there 
is no evidence for significant current star formation in 
the CND of NGC 1068. The starburst activity detected 
by PAH emission likely contributes only 1% o f the to- 
tal infrared luminosity ([Marco fc Brooks! I2OOI and the 
nucle ar starburst is likely 200-300 Myr old (jDavies et al.l 
[2003). 

Oxygen depletion could also be responsible for de- 
creasing the abundance of HCO+, hence increasing the 
HCN/HCO+ ratio. Oxygen underabundance would re- 
duce the formation of CO, leaving more car bon for other 
carbo n ated species (e.g . HCN and HNC. lUsero et al 



TABLE 7 

Z-Spec Determined Line Ratios 
IN NGC1068 



20041) : [Sternberg etaLl ([19941 ): [Shalabiea fc Greenberg 



19961) ). Since HCO+ includes both carbon and oxygen. 



the effects of more available carbon but less available 
oxygen mean no net change in its abundance. There- 
fore, this scenario would predict higher HCN/CO and 
HCN/HCO+ ratio s. This scenario was ruled out by 
[Usero et"al] ([2004D because it predicted a higher ratio of 
HCN/HCO ^ than they f ound: however, our ratio (which 
agrees with [Krips et al.[ ([2008f)') matches bet ter with the 
prediction from the cited [Usero eFall ([2OOI oxygen de- 
pletion model than their cited XDR model. The same is 
true for the HCN/CO abundance ratio. Because their 
oxygen depletion model is based on different physical 
conditions than we find here (nj^j = 1 x 10^ cm~^, T = 20 



Ratio 



Value 



HCN3^2/HNC3^2 
HCO+ 2/HNC3^2 
HCO+_^2/HCN3-,2 

CS6->5/CS5_>4 
CS5_»4/CS4^3 
CS5^4/HNC3^2 

CS5^4/HCO+^2 
CS5^4/HCN3^2 



4.2 ± 0.4 
2.8 ± 0.3 
0.67 ± 0.05 
0.70 ± 0.27 
0.55 ± 0.18 
0.69 ± 0.15 
0.25 ± 0.05 
0.16 ± 0.03 



Note. — Errors do not include 
calibration error (estimated to be 
10%) because they are all measured 
within our Z-Spec band. 



K. IRuffle et"aII(IT998l) ) , we do not comment on the oxy- 
gen depletion scenario based on their models. 

The final possibily is an XDR, which at lower densities 
(such as those inferred from our analysis, median of 10^'^) 
can explain the HCN/HCO"'" ratio. In addition to the 
modeled molecular abundances, the measured line inten- 
sity ratios also lend evidence towards the XDR scenario, 
as shown in the following section. We also found a ratio 
of column densities N(HNC)/N(HCN) < 1 ( 0.3), which 
can be found in XDR environme nts if the total column 
densi ty is lower than 10^^ cm~^ ([Perez-Beaupuits et al.l 
|2003), which we infer through our likelihoods that it is. 

Some other possibilites for the energetics of NGC 1068 
hav e also been c onsidered. The first two have been ruled 
out: [Krips et al. (2008) pointed out that cosmic rays ac- 
celerated in supernova remnants are generally thought to 
increase HCO"*" and decrease HCN abundances, which 
is the opposite of what we observe. IR pumping af- 
fects both HCO+ and HCN, which woul d not explain the 
high HCN/HCO+ ratio. However, Loen en et all (pOOl ) 
found that some XDR results can be explained also by 
a weak PDR with mechanical heating, and the high 
abundance of SiO may ind icate the infiuence of shocks 
([Garcia-Burillo et al.[[20Tot) . 

5.2. Line Ratios: Strong Evidence of XDR Excitation 
of Molecular Gas 

In addition to radiative transfer and non-LTE exci- 
tation modeling, the line ratios of high density tracers 
can themselves be used as diagnostics for XDR envi- 
ronments. Lines within the Z-Spec band benefit from 
our common calibration, making the ratios within our 
band ideal for these diagnostics. Futhermore, even if 
the lines suffer from contamination from the starburst 
ring, if that contamination is of a similar amount for 
the lines we are comparing (i.e. they have similar values 
of fss), then their ratio is unaffected. Some of the in- 
tensity ratios constructed using Z-Spec only are shown 
in Table [Tj T hough our HCO+/ HCN J = 3 ^ 2 ratio 
is higher than [Krips et al.l ([20081) (0.38 ± 0.07) and our 
HCN/HNC slightly lower than [Perez-Beaupuits et al.1 
(j2007.) (6.48 ± 1.95), the ratios cannot be concluded to 
b e incompatible given the er ror bars. 

[Meiierink fc Spaansl ([2005[ ) constructed grids of XDR 
and PDR models. For our molecular lines of interest, 
the grids cover a range of densities (10^ — 10^'^ cm~'^), 
sizes (hence, column densities, 3 x 10^^ — 1 x 10^^ cm^^). 
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and radiation fields. The PDR UV-radiation field is rep- 
resented in multiples of the Habing flux (Gq = 1 corre- 
sponds to 1.6 X 10~^ erg cm~^ s"^) and the XDR X- 
ray field is simply given in erg cm~^ s~^. Intensities 
for various J lines for these molecules are calculated at 
each grid point and contours of some line ratios are pre- 
sented assuming cloud sizes of 1 pc (see Figures 14-17 in 
iMeiierink et al.l (|2007f) ). Assuming this fixed cloud size 
means that the column density varies as a function of 
the density. The line ratios within the Z-Spec band are 
not plotted explicitly in their publication, but the grids 
are available online such that we could make our own 
diagnostic plotfl We constructed 8 line ratios that are 
unique to Z-Spec's band, as presented in Table [3 and 
created contour plots for the XDR and PDR models for 
each ratio. We then constructed a chi-squared statistic 
comparing our ratio to each grid point; the resulting 2- 
sigma confidence intervals (smoothed) for each of the 8 
ratios are presented in Figure [T] 

Of the 8 line ratios we present, 6 are consistent with 
the XDR model at the lower densities inferred by our ra- 
diative transfer modeling (between 10^ and 10^ cm~^). 
We note that the models do not go below 10"* cm~^, 
so for those molecules which favor solutions at this den- 
sity, the parameters should be considered upper limits. 
The remaining 2 ratios are either consistent with XDR 
(HCO+ /HNC J = 3 ^ 2) or PDR (HCO+ /HCN J = 3 2) 
at higher densities. The fact that the two ratios in dis- 
agreement with the others contain HCO"'" may indicate 
something about the physical state of this particular 
molecule; for example, it may be tra cing shock-excited 
material. iGarcia-Burillo et al.l (j2010') have found previ- 
ous evidence of shocks in NGC 1068 using SiO, HNCO, 
and methanol (CH3OH). Our conclusions are similar to 
the aforementioned paper; we find some evidence of XDR 
influence, but also some possible influence of shocks. The 
discrepancy of HCO"*" may also lend support to the oxy- 
gen depletion scenario, because in order to match the 
XDR models at lower densities, both of these ratios 
would need to be significantly higher than we measure. 
If HCN and HNC are being selectively enhanced over 
HCO+, that could drive down these ratios, and hence 
explain why we see lower values than those predicted by 
the XDR models at densities of - lO*-^ cm^^. If HCO+ 
is in fact also tracing gas undergoing a different excita- 
tion mechanism, we again consider the possibility that it 
should not be included in the radiative transfer analysis 
(as in ^4.2p . When we exclude this molecule, the den- 
sity peaks at 10^ ° instead of lO"*'^ cm~'^. Even with the 
slightly larger density, there is no additional support for 
the PDR model, given the arguments already presented 
that remain the same even when HCO"*" is excluded from 
the analysis. However, this is now a second, independent 
line of evidence (in addition to the discussion in 
that suggests that HCO+ may be tracing a different gas 
component than other tracers. 

This analysis adds support to previous studies that 
have demonstrated the XDR nature of t he CND of 
NGC 1068. For example, lUsero et al.l ([200I found sup- 
port for the XDR model based on relative molecular 
abundances derived from LVG simulations (HCN/CO, 

® http:/ /www. strw.lcidcnuniv.nl/~meijerink/grid/ 



CS/C O, HCO/HCO+.CN/HCN). rPerez-Beaupuits et al.l 
(|2009D came to the same conclusion based on line ra- 
tios of HCO+ J = 4 ^- 3 with lower J lines and HCN 
J 4 ^ 3, and by the N(CN)/N(HCN) - 1-4 col- 
umn densi ty ratio. These s tudies are in addition to the 
^arcia -Burillo et al.l (120101 ) conc lusions in the previous 
paragraph. Also. lMalonevI ()1997|) demonstrated that the 
X-ray irradation in the CND is adequate for creating an 
XDR. 

The agreement of the majority of our line intensity ra- 
tios with the XDR scenario at densities consistent with 
our radiative transfer modeling results provides strong 
support for an XDR and also demonstrates the utility of 
Z-Spec's broad bandpass and self-consistent line calibra- 
tion across the band. 

5.3. Black Hole Accretion Rate and CND lifetime 

We used our derived total gas mass (10^ M©) to esti- 
mate the timescale for accretion onto the central black 
hole, based on the accretion rate inferred from the AGN 
luminosity. We use the simple equation 

L^riMc^, (13) 

where ry is the e fficiency. Using an accretion luminosity 
of 3.1 X 10" L^ dRigbv et al.l[200l (based on [OIV] as a 
calibrator) , the accretion rate is 0.2 Mq year~^ assum- 
ing 10% efficiency. Estimating a timescale for accretion 
by McND /M yields ^ 500 Myr. This is only a few times 
the 200-300 Myr estimated age of the nuclear stellar core 
(the innermost 70 pc). By this estimate, the most re- 
cent star formation activity has only taken place for a 
portion of the AGN's depletion timescale. We note that 
the timescale of accretion is proportional to efficiency, 
which is not well known, but most likely we are un- 
derestimating the timescale rather than overestimating, 
considering that the effic iency is likely between 10-100% 
(jSchartmann et al]l2010l ) but not significantly lower than 
10%. There are also uncertainties in the total gas mass, 
le aving it di fficult to precisely compare the timescales. 

iDavies eTal.1 (|2007[ ) found a correlation between the 
age of nuclear starburst and AGN luminosity; of their 
sample of 8 galaxies, NGC 1068 had one of the oldest 
starbursts and also is one of the most luminous AGN, 
accreting at a significant fraction of its Eddington lu- 
minosity. It is interesting that the CND hosts a signif- 
icant amount of dense gas that could be used for star 
formation, yet no significant st ar formation is currently 
happening. IDavies et al.l ()2007D theorize that there is a 
delay between a starburst and AGN activity and that 
the starburst activity can help power the AGN. Though 
the nuclear starburst began 200-300 Myr ago and has 
since effectively ended, its influence is still significant be- 
cause it likely contributed to making NGC 1068 such a 
strong AGN, which now is the primary influence of the 
dense molecular gas as we saw in the XDR analysis of 
the previous section. 

6. CONCLUSIONS AND FUTURE WORK 

We have presented a new broadband spectrum of NGC 
1068 obtained with Z-Spec, complementing the current 
literature with a set of common-calibration measure- 
ments of CO, its isotopes, and other high-density trac- 
ing molecules. After extending the modeling procedure 
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Fig. 7. — 2u confidence intervals for Z-Spec fine ratios. These 
confidence intervals were constructed by comparing Z-Spec's mea- 
sured line ratios (Table [3 with the P DR (top) and XDR (bottom) 
models of lMeiierink fc SpaansI I I2005I '). 



of iWard et al. (I2003D to include multiple species as in 
iNavlor et aiTpOlOO .' we use radiative transfer grids to 
investigate the physical properties of the CND. In our 
preferred model, we correct for a small contribution from 
the starburst ring for IICO+ J = 1 — > but assume that 
all other molecular lines are tracing the CND only. Ra- 
diative transfer analysis cannot substantially constrain 
the temperature of the CND, but we argue there is ev- 
idence for a warm component, so future models should 
reflect higher gas temperatures. Though the tempera- 
ture and density are degenerate, we are able to constrain 
their product (pressure) well, at Log(P [K cm^'^]) = 
6.60 ± 0.24). The column densities and relative abun- 
dances of the molecules, as well as the total mass, are 



also well constrained. We also presented models assum- 
ing greater contributions to the flux by the starburst 
ring, and found that these scenarios are difhcult to re- 
produce with well-constrained models. However, one of 
the models may represent a separate, cooler component 
that could be distinct from warmer gas. Should future 
studies indicate that our preferred model underestimates 
the contribution of the starburst ring, the results would 
change systematically as described in i)4.2l 

In addition to the radiative transfer analysis, we com- 
pared the Z-Spec measured line ratios to those predicted 
by XDR and PDR models. Our line ratios, measured 
using a consistent calibration scheme, match extremely 
well with XDR models at the densities found in our radia- 
tive transfer models. These provide strong independent 
support of the XDR nature of NGC 1068's CND. Both 
the PDR/XDR ratio analysis and the radiative trans- 
fer modeling of HCO+ alone indicate that HCO+ may 
be tracing a different component of gas than the other 
molecules used in this work (HON, HNC, CS). 

Future studies of NGC 1068 by the SPIRE instrument 
on the Herschel Space Observatory will shed further light 
on the conditions of the CND. The Fourier Transform 
Spectrometer should be able to fill out the CO ladder 
from J = 4— )-3 to J = 13— )-12 which will b etter constrain 
the te mperature of the molecular gas (e.g.. lPanuzzo et all 
(|2010f )). Furthermore, the spectroscopic imaging capa- 
bility will offer the opportunity to compare the condi- 
tions of the CND vs. the starburst ring. We predict 
that even some of the higher-J transitions of the high- 
density tracers from the CND will be measureable with 
Herschel's HIFI instrument. ALMA will also be crucial in 
the further study of NGC 1068, especially to resolve the 
separate components within the CND. Higher angular 
resolution is key to understanding the relation between 
the gas and the central source, and also important to 
understanding the clumping and distribution of the gas. 
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